home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / numerical / slatec / zacon.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  10.9 KB  |  273 lines
  1. ;;; Compiled by f2cl version 2.0 beta 2002-05-06
  2. ;;; 
  3. ;;; Options: ((:prune-labels nil) (:auto-save t) (:relaxed-array-decls t)
  4. ;;;           (:coerce-assigns :as-needed) (:array-type ':simple-array)
  5. ;;;           (:array-slicing nil) (:declare-common nil)
  6. ;;;           (:float-format double-float))
  7.  
  8. (in-package "SLATEC")
  9.  
  10.  
  11. (let ((pi_ 3.141592653589793) (zeror 0.0) (coner 1.0))
  12.   (declare (type double-float coner zeror pi_))
  13.   (defun zacon (zr zi fnu kode mr n yr yi nz rl fnul tol elim alim)
  14.     (declare (type (simple-array double-float (*)) yr yi)
  15.              (type f2cl-lib:integer4 kode mr n nz)
  16.              (type double-float zr zi fnu rl fnul tol elim alim))
  17.     (prog ((cyr (make-array 2 :element-type 'double-float))
  18.            (cyi (make-array 2 :element-type 'double-float))
  19.            (cssr (make-array 3 :element-type 'double-float))
  20.            (csrr (make-array 3 :element-type 'double-float))
  21.            (bry (make-array 3 :element-type 'double-float)) (i 0) (inu 0)
  22.            (iuf 0) (kflag 0) (nn 0) (nw 0) (arg 0.0) (ascle 0.0) (as2 0.0)
  23.            (azn 0.0) (bscle 0.0) (cki 0.0) (ckr 0.0) (cpn 0.0) (cscl 0.0)
  24.            (cscr 0.0) (csgni 0.0) (csgnr 0.0) (cspni 0.0) (cspnr 0.0) (csr 0.0)
  25.            (c1i 0.0) (c1m 0.0) (c1r 0.0) (c2i 0.0) (c2r 0.0) (fmr 0.0) (fn 0.0)
  26.            (pti 0.0) (ptr 0.0) (razn 0.0) (rzi 0.0) (rzr 0.0) (sc1i 0.0)
  27.            (sc1r 0.0) (sc2i 0.0) (sc2r 0.0) (sgn 0.0) (spn 0.0) (sti 0.0)
  28.            (str 0.0) (s1i 0.0) (s1r 0.0) (s2i 0.0) (s2r 0.0) (yy 0.0) (zni 0.0)
  29.            (znr 0.0))
  30.       (declare (type (simple-array double-float (2)) cyr cyi)
  31.                (type (simple-array double-float (3)) cssr csrr bry)
  32.                (type double-float znr zni yy s2r s2i s1r s1i str sti spn sgn
  33.                 sc2r sc2i sc1r sc1i rzr rzi razn ptr pti fn fmr c2r c2i c1r c1m
  34.                 c1i csr cspnr cspni csgnr csgni cscr cscl cpn ckr cki bscle azn
  35.                 as2 ascle arg)
  36.                (type f2cl-lib:integer4 nw nn kflag iuf inu i))
  37.       (setf nz 0)
  38.       (setf znr (- zr))
  39.       (setf zni (- zi))
  40.       (setf nn n)
  41.       (multiple-value-bind
  42.           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10
  43.            var-11 var-12)
  44.           (zbinu znr zni fnu kode nn yr yi nw rl fnul tol elim alim)
  45.         (declare
  46.          (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-8 var-9 var-10
  47.           var-11 var-12))
  48.         (setf nw var-7))
  49.       (if (< nw 0) (go label90))
  50.       (setf nn (min (the f2cl-lib:integer4 2) (the f2cl-lib:integer4 n)))
  51.       (multiple-value-bind
  52.           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9 var-10)
  53.           (zbknu znr zni fnu kode nn cyr cyi nw tol elim alim)
  54.         (declare
  55.          (ignore var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-8 var-9 var-10))
  56.         (setf nw var-7))
  57.       (if (/= nw 0) (go label90))
  58.       (setf s1r (f2cl-lib:fref cyr (1) ((1 2))))
  59.       (setf s1i (f2cl-lib:fref cyi (1) ((1 2))))
  60.       (setf fmr (coerce (the f2cl-lib:integer4 mr) 'double-float))
  61.       (setf sgn (coerce (- (f2cl-lib:dsign pi_ fmr)) 'double-float))
  62.       (setf csgnr zeror)
  63.       (setf csgni sgn)
  64.       (if (= kode 1) (go label10))
  65.       (setf yy (- zni))
  66.       (setf cpn (cos yy))
  67.       (setf spn (sin yy))
  68.       (multiple-value-bind
  69.           (var-0 var-1 var-2 var-3 var-4 var-5)
  70.           (zmlt csgnr csgni cpn spn csgnr csgni)
  71.         (declare (ignore var-0 var-1 var-2 var-3))
  72.         (setf csgnr var-4)
  73.         (setf csgni var-5))
  74.      label10
  75.       (setf inu (f2cl-lib:int fnu))
  76.       (setf arg (* (- fnu inu) sgn))
  77.       (setf cpn (cos arg))
  78.       (setf spn (sin arg))
  79.       (setf cspnr cpn)
  80.       (setf cspni spn)
  81.       (if (= (mod inu 2) 0) (go label20))
  82.       (setf cspnr (- cspnr))
  83.       (setf cspni (- cspni))
  84.      label20
  85.       (setf iuf 0)
  86.       (setf c1r s1r)
  87.       (setf c1i s1i)
  88.       (setf c2r (f2cl-lib:fref yr (1) ((1 n))))
  89.       (setf c2i (f2cl-lib:fref yi (1) ((1 n))))
  90.       (setf ascle (/ (* 1000.0 (f2cl-lib:d1mach 1)) tol))
  91.       (if (= kode 1) (go label30))
  92.       (multiple-value-bind
  93.           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9)
  94.           (zs1s2 znr zni c1r c1i c2r c2i nw ascle alim iuf)
  95.         (declare (ignore var-0 var-1 var-7 var-8))
  96.         (setf c1r var-2)
  97.         (setf c1i var-3)
  98.         (setf c2r var-4)
  99.         (setf c2i var-5)
  100.         (setf nw var-6)
  101.         (setf iuf var-9))
  102.       (setf nz (f2cl-lib:int-add nz nw))
  103.       (setf sc1r c1r)
  104.       (setf sc1i c1i)
  105.      label30
  106.       (multiple-value-bind
  107.           (var-0 var-1 var-2 var-3 var-4 var-5)
  108.           (zmlt cspnr cspni c1r c1i str sti)
  109.         (declare (ignore var-0 var-1 var-2 var-3))
  110.         (setf str var-4)
  111.         (setf sti var-5))
  112.       (multiple-value-bind
  113.           (var-0 var-1 var-2 var-3 var-4 var-5)
  114.           (zmlt csgnr csgni c2r c2i ptr pti)
  115.         (declare (ignore var-0 var-1 var-2 var-3))
  116.         (setf ptr var-4)
  117.         (setf pti var-5))
  118.       (f2cl-lib:fset (f2cl-lib:fref yr (1) ((1 n))) (+ str ptr))
  119.       (f2cl-lib:fset (f2cl-lib:fref yi (1) ((1 n))) (+ sti pti))
  120.       (if (= n 1) (go end_label))
  121.       (setf cspnr (- cspnr))
  122.       (setf cspni (- cspni))
  123.       (setf s2r (f2cl-lib:fref cyr (2) ((1 2))))
  124.       (setf s2i (f2cl-lib:fref cyi (2) ((1 2))))
  125.       (setf c1r s2r)
  126.       (setf c1i s2i)
  127.       (setf c2r (f2cl-lib:fref yr (2) ((1 n))))
  128.       (setf c2i (f2cl-lib:fref yi (2) ((1 n))))
  129.       (if (= kode 1) (go label40))
  130.       (multiple-value-bind
  131.           (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9)
  132.           (zs1s2 znr zni c1r c1i c2r c2i nw ascle alim iuf)
  133.         (declare (ignore var-0 var-1 var-7 var-8))
  134.         (setf c1r var-2)
  135.         (setf c1i var-3)
  136.         (setf c2r var-4)
  137.         (setf c2i var-5)
  138.         (setf nw var-6)
  139.         (setf iuf var-9))
  140.       (setf nz (f2cl-lib:int-add nz nw))
  141.       (setf sc2r c1r)
  142.       (setf sc2i c1i)
  143.      label40
  144.       (multiple-value-bind
  145.           (var-0 var-1 var-2 var-3 var-4 var-5)
  146.           (zmlt cspnr cspni c1r c1i str sti)
  147.         (declare (ignore var-0 var-1 var-2 var-3))
  148.         (setf str var-4)
  149.         (setf sti var-5))
  150.       (multiple-value-bind
  151.           (var-0 var-1 var-2 var-3 var-4 var-5)
  152.           (zmlt csgnr csgni c2r c2i ptr pti)
  153.         (declare (ignore var-0 var-1 var-2 var-3))
  154.         (setf ptr var-4)
  155.         (setf pti var-5))
  156.       (f2cl-lib:fset (f2cl-lib:fref yr (2) ((1 n))) (+ str ptr))
  157.       (f2cl-lib:fset (f2cl-lib:fref yi (2) ((1 n))) (+ sti pti))
  158.       (if (= n 2) (go end_label))
  159.       (setf cspnr (- cspnr))
  160.       (setf cspni (- cspni))
  161.       (setf azn (zabs znr zni))
  162.       (setf razn (/ 1.0 azn))
  163.       (setf str (* znr razn))
  164.       (setf sti (* (- zni) razn))
  165.       (setf rzr (* (+ str str) razn))
  166.       (setf rzi (* (+ sti sti) razn))
  167.       (setf fn (+ fnu 1.0))
  168.       (setf ckr (* fn rzr))
  169.       (setf cki (* fn rzi))
  170.       (setf cscl (/ 1.0 tol))
  171.       (setf cscr tol)
  172.       (f2cl-lib:fset (f2cl-lib:fref cssr (1) ((1 3))) cscl)
  173.       (f2cl-lib:fset (f2cl-lib:fref cssr (2) ((1 3))) coner)
  174.       (f2cl-lib:fset (f2cl-lib:fref cssr (3) ((1 3))) cscr)
  175.       (f2cl-lib:fset (f2cl-lib:fref csrr (1) ((1 3))) cscr)
  176.       (f2cl-lib:fset (f2cl-lib:fref csrr (2) ((1 3))) coner)
  177.       (f2cl-lib:fset (f2cl-lib:fref csrr (3) ((1 3))) cscl)
  178.       (f2cl-lib:fset (f2cl-lib:fref bry (1) ((1 3))) ascle)
  179.       (f2cl-lib:fset (f2cl-lib:fref bry (2) ((1 3))) (/ 1.0 ascle))
  180.       (f2cl-lib:fset (f2cl-lib:fref bry (3) ((1 3))) (f2cl-lib:d1mach 2))
  181.       (setf as2 (zabs s2r s2i))
  182.       (setf kflag 2)
  183.       (if (> as2 (f2cl-lib:fref bry (1) ((1 3)))) (go label50))
  184.       (setf kflag 1)
  185.       (go label60)
  186.      label50
  187.       (if (< as2 (f2cl-lib:fref bry (2) ((1 3)))) (go label60))
  188.       (setf kflag 3)
  189.      label60
  190.       (setf bscle (f2cl-lib:fref bry (kflag) ((1 3))))
  191.       (setf s1r (* s1r (f2cl-lib:fref cssr (kflag) ((1 3)))))
  192.       (setf s1i (* s1i (f2cl-lib:fref cssr (kflag) ((1 3)))))
  193.       (setf s2r (* s2r (f2cl-lib:fref cssr (kflag) ((1 3)))))
  194.       (setf s2i (* s2i (f2cl-lib:fref cssr (kflag) ((1 3)))))
  195.       (setf csr (f2cl-lib:fref csrr (kflag) ((1 3))))
  196.       (f2cl-lib:fdo (i 3 (f2cl-lib:int-add i 1))
  197.                     ((> i n) nil)
  198.         (tagbody
  199.           (setf str s2r)
  200.           (setf sti s2i)
  201.           (setf s2r (+ (- (* ckr str) (* cki sti)) s1r))
  202.           (setf s2i (+ (* ckr sti) (* cki str) s1i))
  203.           (setf s1r str)
  204.           (setf s1i sti)
  205.           (setf c1r (* s2r csr))
  206.           (setf c1i (* s2i csr))
  207.           (setf str c1r)
  208.           (setf sti c1i)
  209.           (setf c2r (f2cl-lib:fref yr (i) ((1 n))))
  210.           (setf c2i (f2cl-lib:fref yi (i) ((1 n))))
  211.           (if (= kode 1) (go label70))
  212.           (if (< iuf 0) (go label70))
  213.           (multiple-value-bind
  214.               (var-0 var-1 var-2 var-3 var-4 var-5 var-6 var-7 var-8 var-9)
  215.               (zs1s2 znr zni c1r c1i c2r c2i nw ascle alim iuf)
  216.             (declare (ignore var-0 var-1 var-7 var-8))
  217.             (setf c1r var-2)
  218.             (setf c1i var-3)
  219.             (setf c2r var-4)
  220.             (setf c2i var-5)
  221.             (setf nw var-6)
  222.             (setf iuf var-9))
  223.           (setf nz (f2cl-lib:int-add nz nw))
  224.           (setf sc1r sc2r)
  225.           (setf sc1i sc2i)
  226.           (setf sc2r c1r)
  227.           (setf sc2i c1i)
  228.           (if (/= iuf 3) (go label70))
  229.           (setf iuf -4)
  230.           (setf s1r (* sc1r (f2cl-lib:fref cssr (kflag) ((1 3)))))
  231.           (setf s1i (* sc1i (f2cl-lib:fref cssr (kflag) ((1 3)))))
  232.           (setf s2r (* sc2r (f2cl-lib:fref cssr (kflag) ((1 3)))))
  233.           (setf s2i (* sc2i (f2cl-lib:fref cssr (kflag) ((1 3)))))
  234.           (setf str sc2r)
  235.           (setf sti sc2i)
  236.          label70
  237.           (setf ptr (- (* cspnr c1r) (* cspni c1i)))
  238.           (setf pti (+ (* cspnr c1i) (* cspni c1r)))
  239.           (f2cl-lib:fset (f2cl-lib:fref yr (i) ((1 n)))
  240.                          (- (+ ptr (* csgnr c2r)) (* csgni c2i)))
  241.           (f2cl-lib:fset (f2cl-lib:fref yi (i) ((1 n)))
  242.                          (+ pti (* csgnr c2i) (* csgni c2r)))
  243.           (setf ckr (+ ckr rzr))
  244.           (setf cki (+ cki rzi))
  245.           (setf cspnr (- cspnr))
  246.           (setf cspni (- cspni))
  247.           (if (>= kflag 3) (go label80))
  248.           (setf ptr (coerce (abs c1r) 'double-float))
  249.           (setf pti (coerce (abs c1i) 'double-float))
  250.           (setf c1m (max ptr pti))
  251.           (if (<= c1m bscle) (go label80))
  252.           (setf kflag (f2cl-lib:int-add kflag 1))
  253.           (setf bscle (f2cl-lib:fref bry (kflag) ((1 3))))
  254.           (setf s1r (* s1r csr))
  255.           (setf s1i (* s1i csr))
  256.           (setf s2r str)
  257.           (setf s2i sti)
  258.           (setf s1r (* s1r (f2cl-lib:fref cssr (kflag) ((1 3)))))
  259.           (setf s1i (* s1i (f2cl-lib:fref cssr (kflag) ((1 3)))))
  260.           (setf s2r (* s2r (f2cl-lib:fref cssr (kflag) ((1 3)))))
  261.           (setf s2i (* s2i (f2cl-lib:fref cssr (kflag) ((1 3)))))
  262.           (setf csr (f2cl-lib:fref csrr (kflag) ((1 3))))
  263.          label80))
  264.       (go end_label)
  265.      label90
  266.       (setf nz -1)
  267.       (if (= nw -2) (setf nz -2))
  268.       (go end_label)
  269.      end_label
  270.       (return
  271.        (values nil nil nil nil nil nil nil nil nz nil nil nil nil nil)))))
  272.  
  273.